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ABSTRACT 


A  software  package  for  general  analysis  of  real  time 
series  for  geophysical  data  was  developed.   The  package 
consists  of  FORTRAN  IV  callable  subprograms  that  are  as- 
signed different  tasks  such  as  input,  output  (either 
printed  or  graphical) ,  initial  data  analysis  and  filter- 
ing ,  spurious  data  rejection,  FFT,  and  spectral  analysis. 
A  geophysical  problem  involving  detection  of  unusual  deep 
ocean  pressures  was  analyzed. 
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I.   INTRODUCTION 

The  types  of  real  time  series  considered  here  are  geo- 
physical processes.   Usually,  time  is  the  independent  vari- 
able and  such  quantities  as  velocity  and  direction  of  flow, 
temperature,  and  pressure  are  the  dependent  variables.   In 
the  analysis  and  design  of  geophysical  experiments,  it  is 
useful  to  have  at  one's  disposal  a  set  of  software  sub- 
routines to  aid  in  this  analysis.   Thus,  a  software  pack- 
age for  general  analysis  of  real  time  series  geophysical 
data  was  developed.   The  package  consists  of  FORTRAN  IV 
callable  subprograms  that  are  assigned  different  tasks  such 
as  input,  output  (either  printed  or  graphical),  initial 
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spectral  analysis,  and  so  forth.  A  geophysical  problem 
involving  detection  of  unusual  deep  ocean  pressures  was 
analyzed . 


II.   A  GEOPHYSICAL  PROBLEM 

A.  STATEMENT  OF  THE  PROBLEM 

Throughout  nature  there  exist  many  phenomenon  that 
occur  only  as  "unusual  events."   In  the  realm  of  oceano- 
graphy one  such  unusual  event  is  called  a  tsunami.    The 
word  tsunami  comes  from  the  Japanesse  and  means  storm  wave. 
It  is  known  that  the  presence  of  a  tsunami  is  indicated  by 
a  definite  change  in  water  pressure.   Other  information 
concerning  the  velocity  and  direction  of  the  current,  water 
temperature,  salinity,  and  so  forth,  would  also  be  useful 
in  the  analysis  and  understanding  of  these  events.   However, 
the  infrequent   rate  at  which  tsunamis  occur  makes  their 
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these  events,  methods  were  devised  to  detect  and  record 
them.   The  problem  then  is  to  logically  design  a  detector 
that  will  automatically  sense  the  presence  of  a  tsunami 
and  cause  all  the  pertinent  data  to  be  recorded. 

B.  A  MODEL 

In  order  to  determine  the  presence  of  an  unusual  event, 
a  detector  which  receives  digitized  pressure  data  was  de- 
signed [1]  .   Throughout  the  following  discussion  refer  to 
Figure  1  for  the  logical  design  of  the  detector.   The  de- 
tector can  be  divided  into  six  logically  separate  sections. 

Data  Input.   Digital  pressure  sensors  of  varying 
degrees  of  sophistication  and  reliability  are  presently 
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available  on  the  commercial  market.   Digital  sensors 
generally  count  the  number  of  cycles  of  a  pressure  sensi- 
tive oscillator  over  a  given  time  interval.   This  summa- 
tion in  effect  integrates  the  pressure  signal.   The 
integration  process  serves  as  a  low  pass  filter.   That 
is,  high  frequency  pressure  variations  relative  to  the 
sampling  interval  tend  to  average  out,  thus  contributing 
equally  to  each  sample.   This  filtering  of  high  frequency 
variations  might  prove  undesirable  in  some  applications. 
For  this  particular  problem,  however,  almost  all  of  the 
high  frequency  variations  are  naturally  filtered  out  be- 
cause of  attenuation  of  pressure  variations  due  to  ocean 
depth. 

It  was  determined  thai  the  frequency  band  in  which 
tsunamis  occur  extends  from  \    to  30  cycles  per  hour. 
Nyquist  frequency  theory  indicates  that  in  order  to  detect 
a  frequency  of  N  cph ,  a  sampling  rate  of  at  least  2N  cph 
is  required.   Hence  the  sampling  interval  was  determined 
to  be  one  sample  per  minute. 

Input  Buffer.   An  input  buffer  receives  the  most  re- 
cent pressure  reading  in  digitized  form.   The  actual  hard- 
ware implementation  of  the  input  buffer  is  not  part  of  the 
problem.   Conceivably  the  input  buffer  could  be  anything 
from  an  acoustic  delay  line  to  a  magnetic  core. 

The  buffer  is  large  enough  to  retain  the  most  recent 
121  consecutive  samples.   Hence,  a  two  hour  history  of 
pressure  data  is  available  for  analysis  at  any  one  time. 


Figure  2(a)  is  a  plot  of  typical  pressure  data  over  an  80 
hour  period. 

First  Differencing.   The  presence  of  a  tsunami  could 
be  indicated  by  an  abnormal  change  in  pressure  from  one 
point  in  time  to  another.   The  phrase  "change  in  pressure" 
indicates  that  the  derivative  should  be  computed.   As  a 
close  approximation  to  the  derivative,  the  positive  first 
difference  was  determined. 

Incorporated  within  the  first  differencer  was  a  spur- 
ious data  gate.   Here  each  first  difference  was  compared 
against  a  preset  threshold  value.   If  the  threshold  value 
was  exceeded,  then  a  sensor  error  (spurious  error)  was  the 
probable  cause  since  deep-ocean  pressure  changes  are  rela- 
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difference  might  falsely  indicate  the  presence  of  a  tsunami, 
an  excessive  first  difference  was  set  equal  to  the  previous 
first  difference.   The  threshold  value  was  determined  after 
analysis  of  test  data.   First  differencing,  as  computing 
the  first  derivative,  removes  the  mean  pressure  value  and 
suppresses  drift.   Figure  2(b)  illustrates  the  effect  of 
first  differencing  on  the  data  of  Figure  2(a). 

Tide  Filter.   The  purpose  of  the  tide  filter  is  to 
remove  the  effects  of  the  diurnal  and  semi-diurnal  tides. 
This  was  accomplished  by  forming  a  weighted  sum  from  the 
sample  data  at  h   hour  intervals.   The  coefficients  were 
determined  from  analysis  of  the  tidal  spectrum.   The  tide 
filter  output  then  is: 
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D«.  =  dP+  -  3.920392-dP,.  xn    +    5 . 841800 -dP+  ,n 
t  t  t-30       •       t-60 


-  3. 920392- dPt.9Q  +  dPt.120 


where  dP   is  the  first  difference  at  time  t  [2] . 

Figure  2(c)  shows  the  effect  of  the  tide  filter  on 
the  first  differenced  data  of  Figure  2(b).   Note  how  the 
tide  filter  extended  the  effects  of  a  spurious  data  value 
when  there  was  no  spurious  data  gating. 

Rectification.   As  Figure  2(c)  shows,  the  output  of 
the  tide  filter  presents  an  erratic  picture  of  the  changes 
in  pressure. ■  In  order  to  observe  any  trends  in  this  data, 
an  integral  over  some  fixed  period  of  time  needed  to  be 
computed.   Any  integration  of  the  tide  filtered  data 
would  only  produce  a  result  close  to  zero.   This  problem 
was  overcome  by  taking  the  absolute  value  of  the  data  in 
two  ways;   that  is,  the  negative  data  values  were  made 
positive  so  that  the  data  was  no  longer  distributed  about 
zero.   The  most  straightforward  method  was  to  simply  ignore 
the  sign  (i.e.,  rectification)  of  each  number  outputted 
from  the  tide  filter.   The  second  method  was  to  square  each 
value.   Figure  2(d)  shows  the  result  of  squaring  the  tide 
filter  output. 

Integration  and  Detection.   The  rectified  or  squared 
tide  filtered  data  was  then  integrated.   The  method  of 
integration  was  based  on  the  solution  of  a  differential 
equation.   The  general  form  of  the  integrating  function  is: 
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y.+-,  =  e  a  (yt+xt)  >   t  =  1,2,3,...  minutes 

and  where  y,  =  0 , 

a  =  a  time  constant  associated  with  the  inter- 
val of  integration, 
x   =  the  data  to  be  integrated  (rectified  or 
squared  output  of  tide  filter) . 
The  value  of  a  was  determined  experimentally  via  the  simu- 
lation (a  =  1/60).   Figure  2(e)  is  a  plot  of  the  integrated 
tide  filtered  pressure  data  after  the  squaring  process. 
Dependent  on  the  value  of  a,  a  critical  threshold 
value  had  to  be  chosen  such  that  whenever  the  integrated 
signal  exceeded  this  value,  a  switch  was  closed  causing 
the  contents  of  the  buffer  to  be  recorded  on  a  tape  or 
otherwise  be  permanently  recorded.   At  the  end  of  a  two 
hour  period  of  normal  input  the  buffer  would  once  again 
be  recorded.   At  this  point  the  value  of  the  integrated 
signal  would  be  tested  again.   If  the  integrated  signal 
was  still  greater  than  the  threshold  value,  then  the  buf- 
fer would  again  be  recorded  after  two  hours;  otherwise 
after  two  hours  the  testing  is  resumed  as  before.   This 
technique  guarantees  that  a  minimum  four  hour  record  will 
be  available  for  each  event  detected. 
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FIGURE  2. 


(a)  Pressure,  (b)  First  Differenced  Pressure  (spur 
ious  Data  Gate  Off),  (c)  Tide  Filtered  Pressure, 
(d)   Squared  Tide  Filtered  Data,  (e)  Integrated 
Signal  Versus  Time  in  Hours. 
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III.   A  SOFTWARE  SIMULATION 

A.    GENERAL  APPROACH 

To  be  generally  useful  a  software  package  must  be  able 
to  simulate  a  real  system  on  the  level  of  logical  components. 
That  is,  the  basic  information  extracting  procedures  of  the 
system  must  be  available  for  graphical  or  statistical  anal- 
ysis during  the  simulation.   Accordingly,  the  following  set 
of  eighteen  FORTRAN  callable  subprograms  was  developed. 
The  subprograms  and  their  associated  formal  parameter  are 
listed  alphabetically  in  Table  I.   The  subprogram  UNIFRM 
is  a  function  subprogram  and  XYAXIS  is  an  entry  in  DRAW1 . 
Otherwise  all  the  subprograms  listed  in  Table  I  are  sub- 
routines.  The  various  subprograms  are  explained  in  the 
following  sections  according  to  their  use  in  the  unusual 
event  problem.   Each  subprogram  is  listed  in  its  entirety 
in  the  Appendix. 

Each  subroutine  uses  the  variable  FILE  as  its  first 
formal  parameter.   The  meaning  and  use  of  FILE  will  be 
given  once  here  to  avoid  duplication  in  the  sections  which 
follow. 

FILE  is  a  real-valued  variable  identifying  the  current 
data  file.  This  number  is  for  user  identification  purposes 
only.  FILE  is  usually  set  equal  to  1.0  or  the  DO-loop  in- 
dex if  the  simulation  is  being  performed  iteratively.  When 
a  subroutine  is  called,  a  message  is  printed  indicating 
that  the  subroutine  has  completed  its  task.  The  message  is 
of  the  form: 
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SUBROUTINE 

FORMAL  PARAMETERS 

PAGE 

ANT IDE 

(FILE, 
(FILE, 

Y, 
Y, 

TIDE,  L) 
L) 

26 

CARDS 

16 

DRAW1 

(FILE, 
(FILE, 

Y. 

A. 

L,  YMAX,  YMIN,  LOG,  XL, 
B,  N,  ISN) 

XI)   29 

FFREAL 

23 

FFT 

(FILE, 

A; 

B,  N,  ISN) 

23 

FIRSTD 

(FILE, 

Y; 

DIFF,  L,  THRESH) 

25 

HARMON 

(FILE, 

Y 

,  YH,  L,  HI,  H2) 

24 

INTGRL 

(FILE 

x, 

Y,  L,  FTC) 

27 

OUTPUT 

(FILE, 

Y 

.  L) 

28 

PACK 

(FILE, 
ISEED 

Y 
) 

,  L,  Y2,  L2,  IDEN,  RAMP, 

19 

POWER 

(FILE, 

Y. 

L,  XL,  XI) 

25 

RECSQR 

(FILE,  Y 

(FILE,  Y 

(FILE,  Y 
YSTDEV) 

(FILE,  Y 

(ISEED) 

(XAL,  YA1 

v  9    t     t  n  c  "\ 

f                   l        w   J       i-J   ?       ±    1\C>    J 

,  L,  Y2,  L2) 

,  L,  YMEAN,  YSTDEV,  L2) 

,  L,  YMAX,  YMIN,  YMEAN, 

,  L) 

0 

26 

RESET 

18 

SCAN 

.  17 

STATS 

17 

TAPE 

16 

UNIFRM 

22 

XYAXIS 

29 

TABLE  I. 


Subroutine  Names,  Formal  Parameters,  and  Page  References 
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FILE  NO.  XX. X  "message". 

The  digits  XX. X  indicate  the  current  value  of  FILE.   Execu- 
tion of  subroutine  DRAW1  causes  FILE  to  be  incremented  by 
0.1. 

B.    DATA  INPUT 

CALL  CARDS     (FILE,  Y,  L) 

CALL  TAPE      (FILE,  Y,  L) 

Data  input  is  accomplished  by  calling  one  of  the  above 
subroutines.   The  meaning  of  each  of  the  formal  parameters 

« 

is  : 

FILE    -  see  Section  III. A. 

Y      -  A  real,  singly  dimensioned  array  into  \\rhich 
the  data  is  to  be  read. 

L       -  An  integer- valued  variable  equal  to  the  size 

of  array  Y.   In  subroutine  CARDS,  L  must  equal 
the  number  of  data  elements  to  be  read-in.   In 
subroutine  TAPE,  L  may  equal  the  number  of 
data  elements  to  be  read-in,  or  it  may  equal 
the  maximum  size  of  Y  as  defined  in  the  DIMEN- 
SION statement  of  the  main  program.   If  an 
end-of-file  mark  is  encountered  on  the  tape 
before  the  L  data  elements  have  been  read,  the 
subroutine  halts  and  returns  the  number  of  ele 
ments  actually  processed  as  the  value  of  L. 

These  subroutines  should  be  reloaded  with  the  proper 
FORMAT  statements  inserted  before  the  simulation  is  run. 
Currently  the  FORMAT  is  (10F7 . 0  ,F10 . 0)  . 
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C.    INITIAL  SCAN  OF  DATA 

CALL  STATS    (FILE,  Y,  L,  YMAX ,  YMIN,  YMEAN,  YSTDEV) 
CALL  SCAN     (FILE,  Y,  L,  YMEAN,  YSTDEV,  L2) 
It  is  sometimes  the  case  when  geophysical  data  is  being 
recorded  that  the  initial  start  up  of  the  sensing  device 
intorudces  spurious  data  which  is  recorded  as  valid  data. 
Used  in  series ,  subroutines  STATS  and  SCAN  can  determine 
the  extent  of  any  start  up  noise  present  in  the  data. 
The  meaning  of  each  of  the  formal  parameters  is: 
FILE    -  see  Section  III. A. 

Y      -  A  real-valued,  singly  dimensioned  array  con- 
taining the  data  to  be  analyzed.   Neither 

subroutine  alters  the  contents  of  array  Y  in 

any  way . 
L      -  An  integer-valued  variable  equal  to  the  size 

of  array  Y. 
YMAX    -  A  real-valued  number  returned  by  subroutine 

STATS  equal  to  the  maximum  value  encountered 

in  array  Y. 
YMIN    -  A  real-valued  number  returned  by  subroutine 

STATS  equal  to  the  minimum  value  encountered 

in  array  Y. 
YMEAN   -  A  real-valued  number  returned  by  subroutine 

STATS  equal  to  the  statistical  mean  of  the 

data  in  array  Y. 
YSTDEV  -  A  real-valued  number  returned  by  subroutine 

STATS  equal  to  the  standard  deviation  of  the 

data  in  array  Y. 
17 


L2      -  An  integer-valued  number  returned  from  sub- 
routine SCAN  equal  to  the  number  of  noise 
dependent  elements  in  array  Y  as  determined 
from  the  following  procedure. 

1)  set  1=1,  set  ITEST=0. 

2)  If  ABS(Y(I)-YMEAN)>2*STDEV  then  set 

ITEST=0,  set  L2=I ; 
otherwise,  set  ITEST=ITEST+1 . 

3)  If  ITEST  >  100  then  return;  otherwise 

set  1=1+1,  go  to  step  2. 

D.    ADJUSTING  THE  SAMPLING  INTERVAL 

CALL  RESET    (FILE,  Y,  L,  Y2 ,  L2) 

CALL  PACK     (FILE,  Y,  L,  Y2 ,  L2,  IDEN,  RAMP,  ISEED) 

UNIFRM       (ISEED) 

The  sample  size  of  data  to  be  analyzed  can  be  adjusted 
in  two  basic  ways.   First,  the  sampling  interval  may  be 
shifted  forward  by  any  specified  amount.   For  example,  if 
preliminary  analysis  of  the  data  indicates  that  the  initial 
data  elements  are  in  error,  then  the  data  file  can  be  re- 
set so  that  the  record  begins  later  in  the  sequence,  there- 
by skipping  the  bad  data.   This  is  accomplished  by  a  call 
on  subroutine  RESET.   The  formal  parameters  of  RESET  are: 

FILE    -  see  Section  III. A. 

Y      -  A  real-valued,  singly  dimensioned  array  con- 
taining the  data  to  be  reset. 

L      -  An  integer-valued  variable  equal  to  the  size 
of  array  Y. 
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Y2      -  A  real-valued,  singly  dimensioned  array  into 
which  the  reset  elements  of  array  Y  are  to  be 
loaded.   Arrays  Y  and  Y2  may  be  the  same 
physical  array  in  which  case  the  skipped  data 
is  overlayed  and  lost  to  further  processing. 
L2      -  An  integer-valued  variable.   On  the  initial 
call  to  RESET,  L2  equals  the  number  of  data 
elements  in  array  Y  which  are  to  be  skipped. 
On  return,  L2  equals  the  size  of  array  Y2. 
The  second  method  of  adjusting  the  sampling  interval 
involves  increasing  the  density  of  the  data  points.   This 
is  accomplished  by  a  call  on  subroutine  PACK.   The  formal 
parameters  of  PACK  are: 

FT  f  P       _    c  e\  r*   C  rt  r  -f  i  Ap  A 

Y       -  A  real-valued,  singly  dimensioned  array  con- 
taining the  data  to  be  packed. 

L      -  An  integer-valued  variable  equal  to  the  size 
of  array  Y. 

Y2      -  A  real-valued,  singly  dimensioned  array  into 
which  the  expanded  (packed)  data  of  array  Y 
is  to  be  placed.   Arrays  Y  and  Y2  may  be  the 
same  physical  array. 

L2      -  An  integer- valued  variable.   On  the  initial 
call  to  the  subroutine,  L2  must  equal  the 
maximum  size  of  array  Y2  as  declared  in  the 
dimension  statement  of  the  main  program.   On 
return,  L2  will  equal  the  size  of  array  Y2 
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as  declared  in  the  dimension  statement  of 
the  main  program.   On  return,  L2  will  equal 
the  size  of  array  Y2  if  the  subroutine  com- 
pleted successfully.   The  size  of  array  Y2 
can  be  determined  from  the  formula  (L-l)x 
IDEN+1. 

IDEN    -  An  integer-valued  variable  specifying  the 

multiple  of  the  current  density  that  is  de- 
sired.  The  value  of  IDEN  in  part  determines 
the  value  of  L2  that  is  returned  to  the  main 
program.   IDEN  can  also  be  interpreted  to 
mean  that  IDEN-1  points  are  embedded  between 
each  of  the  original  points  in  array  Y. 

PwAM^    -  \   v>r> al  - v alue d  ,rari?ble  e^ual  to  the  maximum 
variability  allowed  when  the  embedded  points 
are  being  computed.   If  RAMP  =  0.0  then  the 
embedded  points  are  determined  by  linear  in- 
terpolation between  successive  Y  values; 
otherwise,  the  value  of  RAMP  is  multiplied 
by  a  uniformly  random  number  between  -%,   and 
^,  and  this  value  is  added  to  the  interpolated 
value . 

ISEED   -  An  integer-valued  variable  that  serves  as  the 
random  number  seed  for  a  random  number  gener- 
ator.  ISEED  should  be  a  six  digit  integer 
equal  to  ±3(mod  8). 

Conceptually  subroutine  PACK  embeds  the  data  as  shown 
in  Figure  3. 
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FIGURE  3. 


Conceptual  Operation  of  Subroutine  PACK  where  RAMP  >  0 
and  I  DEN  =  4. 
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Function  subprogram  UNIFRM  returns  a  uniformly  random 
number  between  -\   and  ^.   The  use  of  UNIFRM  is  not  re- 
stricted to  subroutine  PACK. 

The  original  data  used  in  the  simulation  was  based  on 
a  four  minute  sampling  interval.   The  detector,  however, 
was  designed  to  operate  on  one  minute  data.   Hence,  to  ob- 
tain the  proper  sampling  interval  the  following  call  was 
made  on  subroutine  PACK. 

CALL  PACK  (FILE,  Y,  L,  Y,  LMAX ,  4,  0.0,  103941) 

LMAX  equals  the  maximum*  size  of  array  Y  as  dimensioned  in 
the  main  program.   Note  RAMP  equals  zero. 

E.    FOURIER  ANALYSIS 

TAT T      pFRFAT  fPTLE        ARM        TSN^ 

CALL  FFT       (FILE,  A,  B,  N,  ISN) 

CALL  HARMON    (FILE,  Y,  YH ,  L,  HI,  H2) 

CALL  POWER     (FILE,  Y,  L) 

Three  areas  of  Fourier  analysis  are  provided  for  in 
the  software  package.   Subroutine  FFREAL  returns  the  Fourier 
cosine  and  sine  coefficients  of  the  data  passed  to  it.   Sub- 
routine HARMON  returns  the  approximation  to  the  data  passed 
to  it  on  the  basis  of  the  HI  through  H2  harmonics.   Sub- 
routine POWER  returns  the  power  spectrum  of  the  data.   Sub- 
routine HARMON  and  POWER  both  call  on  FFREAL,  and  FFREAL  in 
turn  calls  on  subroutine  FFT.   Subroutines  FFREAL  and  FFT 
are  variations  of  Singleton's  [3]  algorithms  for  mixed  radix 
fast  Fourier  transform. 
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The  formal  parameters  of  subroutine  FFT  are: 

FILE    -  see  Section  III. A. 

A,  B    -  Real-valued,  singly  dimensioned  arrays.   Ar- 
rays A  and  B  originally  hold  the  real  and 
imaginary  parts  of  the  data  and  return  the 
real  and  imaginary  Fourier  coefficients. 

N      -  An  integer-valued  variable  equal  to  the  num- 
ber of  real  (or  imaginary)  data  elements. 
In  the  context  of  this  paper,  N=L/2. 

ISN     -  An  integer- valued  variable.   The  magnitude  of 
ISN  determined  the  step  size  of  data  in  ar- 
rays A  and  B.   The  sign  of  ISN  determines 
whether  the  complex  or  inverse  transform  will 
be  performed. 

The  formal  parameters  for  subroutine  FFREAL  are: 

FILE    -  see  Section  III. A. 

A,  B,   -  Real-valued,  singly  dimensioned  arrays.   The 
original  data  values  are  alternately  stored 
in  arrays  A  and  B.   On  return,  arrays  A  and 
B  hold  the  cosine  and  sine  coefficients  re- 
spectively, of  the  transform. 

N      -  An  integer- valued  variable  equal  to  the  num- 
ber of  elements  in  array  A.   In  the  context 
of  this  paper,  N=L/2. 

ISN     -  An  integer- valued  variable.   The  magnitude  of 
ISN  determines  the  step  size  within  arrays  A 
and  B.   If  the  magnitude  of  ISN  is  1,  then 
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each  consecutive  element  of  the  array  is  used 
in  the  computations.   If  the  magnitude  is  2, 
every  other  element  is  used.   Hence,  if  the 
data  is  contained  within  a  single  array,  say 
X,  rather  than  two  separate  ones,  the  call 
to  FFREAL  would  be 

CALL  FFREAL  (X,  X(2) ,  N,  2). 

If  the  sign  of  ISN  is  negative,  then  the  in- 
verse transform  is  computed. 

The  formal  parameters  of  subroutine  HARMON  are: 

FILE    -  see  Section  III. A. 

Y      -  A  real-valued,  singly  dimensioned  array  orig- 

XUdxx^       ^OiH-aillliig       Cii^       v-t  a.  i,  a.       t~  <->      u  L,      aiia-i.)   *_ »_  ti  . 

On  return  array  Y  contains  the  Fourier  co- 
sine and  sine  coefficients  as  returned  from 
FFREAL. 

YH     -  A  real-valued,  singly  dimensioned  array  which 
on  return  holds  the  approximation  based  on 
the  HI  and  H2  harmonic  of  the  data  originally 
contained  in  array  Y.   This  process  involves 
a  double  summation  and  could  be  time  consum- 
ing depending  on  the  values  of  HI  and  H2. 

L      -  An  integer-valued  variable  equal  to  the  size 
of  array  Y  (or  array  YH) . 

HI,  H2  -  Positive  integer- valued  variables  indicating 
which  harmonics,  HI  through  H2,  inclusive,  of 
array  Y  are  to  be  isolated. 
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The  formal  parameters  of  subroutine  POWER  are: 
FILE    -  see  Section  III. A. 

Y  -  A  real-valued,  singly  dimensioned  array  con- 

taining the  data  which  is  to  have  its  power 
spectrum  computed.   On  return,  the  power 
spectrum  coefficients  are  in  the  first  N=L/2 
elements  of  array  Y. 
L       -  An  integer-valued  variable  equal  to  the  size 
of  array  Y. 

F.    FIRST  DIFFERENCING 

CALL  FIRST.D    (FILE,  Y,  DIFF,  L,  THRESH) 

Subroutine  FIRSTD  computes  the  positive  first  dif- 
ference of  the  data  in  array  Y,  and  returns  the  array  of 
first  differences  in  array  DIFF.   The  first  and  last  ele- 
ments of  DIFF  are  zero.   Subroutine  FIRSTD  is  so  coded  that 
arrays  Y  and  DIFF  may  be  the  same  physical  array.   By  set- 
ting THRESH  to  any  computed  value,  any  first  difference 
greater  than  THRESH  will  cause  that  first  difference  to  be 
set  equal  to  the  previous  first  difference. 

The  formal  parameters  of  subroutine  FIRSTD  are: 

FILE    -  see  Section  III. A. 

Y  -  A  real-valued,  singly  dimensioned  array  con- 

taining the  data  to  be  first  differenced. 
DIFF    -  A  real-valued,  singly  dimensioned  array  which 
returns  the  positive  first  differences  of  the 
data  in  array  Y. 
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L      -  An  integer-valued  variable  equal  to  the  size 
of  array  Y  (or  array  DIFF) . 

THRESH  -  A  real-valued  variable  against  which  each 

first  difference  is  compared.   If  the  current 
first  difference  is  greater  than  THRESH 
(threshold  value) ,  then  the  first  difference 
is  set  equal  to  the  previous  first  difference. 

G.    TIDE  FILTER 

CALL  ANTIDE    (FILE,  Y,  TIDE,  L) 

Geophysical  events  usually  occur  in  the  domain  of  other 
less  interesting  phenomenon.   In  the  case  of  the  unusual 
event  detector,  one  such  phenomenon  is  the  effect  of  the 
tides,  both  diurnal  and  semi-diurnal.   Subroutine  ANTIDE 
(anti-tide)  filters  out  the  combined  effect  of  these  tides. 
The  formal  parameters  of  subroutine  ANTIDE  are: 
FILE    -  see  Section  III. A. 

Y  -  A  real-valued,  singly  dimensioned  array  con- 
taining the  data  to  be  tide  filtered. 
TIDE  -  A  real-valued,  singly  dimensioned  array  into 
which  the  tide  filtered  data  is  placed.  Ar- 
rays Y  and  TIDE  may  be  the  same  physical  ar- 
ray. 
L  -  An  integer- valued  variable  equal  to  the  size 
of  array  Y  (or  array  TIDE) . 

H.    RECTIFYING  THE  SIGNAL 

CALL  RECSQR    (FILE,  Y,  Y2 ,  L,  IRS) 

Two  methods  of  rectifying  the  data  are  provided.   The 
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first  involves  taking  the  absolute  value  of  each  data  value; 
while  the  second  squares  each  data  value. 

The  formal  parameters  of  subroutine  RECSQR  (rectify- 
square)  are  : 

FILE    -  see  Section  III. A. 

Y      -  A  real-valued,  singly  dimensioned  array  con- 
taining the  data  to  be  rectified  or  squared. 

Y2      -  A  real-valued,  singly  dimensioned  array  into 
which  the  rectified  or  squared  data  is  to  be 
placed.   Arrays  Y  and  Y2  may  be  the  same 
physical  array. 

L       -  An  integer-valued  variable  equal  to  the  size 
of  array  Y  (or  array  Y2) . 

IRS     -  An  in  Lege i- valued  variable  equal  to  1  or  2. 
If  IRS  =  1,  then  the  data  will  be  rectified 
by  computing  the  absolute  value  of  each  ele- 
ment.  If  IRS  =  2,  then  the  data  will  be 
rectified  by  computing  the  square  of  each 
element . 

I.    INTEGRATION 

CALL  INTGRL    (FILE,  X,  Y,  L,  FTC) 

Integration  is  performed  by  multiplying  the  sum  of 
previous  data  elements  by  a  fading  time  constant  (FTC) ; 
that  is  some  positive  number  less  than  unity.   The  form  of 
the  integrating  equation  is: 

Y(I)  =  EXP(-FTC)  *  (  Y(I-l)  +  X(I)  ) 
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where  FTC  is  the  fading  time  constant,  and  I  is  the  index 
of  the  current  variable. 

The  formal  parameters  of  subroutine  INTGRL  are: 

FILE    -  see  Section  III. A. 

X      -  A  real-valued,  singly  dimensioned  array  con- 
taining the  data  to  be  integrated. 

Y  -  A  real-valued,  singly  dimensioned  array  into 

which  the  integrated  data  is  placed. 
L       -  An  integer-valued  variable  equal  to  the  size 

of  array  X  (or  array  Y) . 
FTC     -  A  real-valued  variable.   A  sample  value  of 

FTC  could  be  1/60. 

J.    OUTPUT 

CALL  OUTPUT    (FILE,  Y,  L) 

CALL  DRAW1     (FILE,  Y,  L,  YMAX,  YMIN,  LOG,  XL,  XI) 

CALL  XYAXIS    (XAL,  YAL) 

Two  means  of  output  are  provided.   Subroutine  OUTPUT 
prints  out  the  first  L  elements  of  array  Y.   The  page  for- 
mat of  OUTPUT  is  200  numbers  with  index  numbers  per  page. 

The  formal  parameters  of  subroutine  OUTPUT  are: 

FILE    -  see  Section  III. A. 

Y  -  A  real-valued,  singly  dimensioned  array  con- 

taining the  data  to  be  printed. 
L       -  An  integer-valued  variable  equal  to  the  num- 
ber of  elements  of  array  Y  that  are  to  be 
printed.   L  may  equal  the  full  size  of  array 
Y  in  which  case  a  large  quantity  of  paper  is 
usually  generated. 
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Subroutine  DRAW1  plots  the  first  L  elements  of  array 
Y.   The  formal  parameters  of  subroutine  DRAW1  are: 

FILE    -  see  Section  III. A. 

Y      -  A  real-valued,  singly  dimensioned  array  con- 
taining the  data  to  be  plotted.   Subroutine 
DRAW1  does  not  alter  the  contents  of  array  Y. 

L      -  An  integer-valued  variable  equal  to  the  num- 
ber of  data  points  to  be  plotted. 

YMAX,   -  Real-valued  variables  equal  to  the  maximum 

YMIN     and  minimum  values  to  be  plotted.   If  YMAX  < 
YMIN,   then  the  plot  is  autoscaled. 

LOG     -  An  integer -valued  variable  equal  to  0  or  1 . 

If  LOG  =  1,  then  the  data  in  array  Y  is  plot- 

i,V_  U.    VJll    J-  \_/  g  .  .     JV^CXXO.       NjcilV-.Iwi.JV_.,     A.  J-    l_IW^  \J 

the  data  is  plotted  in  the  normal  fashion. 
XL,  XI  -  Real-valued,  singly  dimensioned  arrays  con- 
taining literal  constants.   The  maximum  length 
string  allowable  is  40  characters.   XL  is  the 
plot  label  which  may  be  blank,  and  XI  is  the 
user's  identification  label. 
The  lengths  of  the  plot  axes  are  initialized  to  20 
inches  for  the  x-axis  and  8  inches  for  the  y-axis.   However, 
these  dimensions  may  be  changed  by  calling  XYAXIS.   The  for- 
mal parameters  of  XYAXIS  are: 

XAL     -  A  real-valued  variable  equal  to  the  new  length 
of  the  x-axis  in  inches. 
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YAL     -  A  real-valued  variable  equal  to  the  new  length 
of  the  y-axis  in  inches.   The  value  of  YAL 
should  be  less  than  9  inches. 
Subroutine  DRAW1  (and  XYAXIS)  make  use  of  the  standard 
plot  package  for  the  CALCOMP  765  plotter  available  on  the 
IBM  360  at  the  Naval  Postgraduate  School. 
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IV.   CONCLUSIONS 

In  the  detector,  the  question  arose  as  to  what  was  the 
best  value  of  the  fading  time  constant  used  in  the  integra- 
tion scheme.   That  is,  given  the  form  of  the  integrating 
equation  and  that  integration  is  naturally  a  low  pass  fil- 
ter, what  value  of  the  fading  time  constant  would  produce 
the  most  favorable  signal-to-noise  ratio?   Three  values  of 
the  time  constant  were  chosen,  and  the  integration  was  per- 
formed three  times  on  the  squared  output  of  the  tide  filter 
with  a  different  time  constant  in  effect  each  time.   The 
value  of  the  constant  used  to  produce  Figure  2  was  1/60. 
The  significance  of  this  value  is  that  because  of  the  struc- 
ture of  the  intear?-tina  equation   the  contribution  of  a 
single  data  value  after  one  hour's  integration  would  have 
been  reduced  to  one  third  of  its  original  effect.   This 
value  of  the  time  constant  was  used.   The  other  two  values 
of  the  time  constant  considered  were  1/30  and  1/120. 
Table  II  shows  the  results  of  two  signal-to-noise  ratio 
computations  using  the  three  constants.   The  two  hour 
"threshold"  level  to  noise  ratio  was  computed  by  finding  a 
level  such  that  the  signal  exceeded  that  level  for  exactly 
two  hours . 

The  integrations  were  performed  a  second  time  but  the 
rectified  rather  than  the  squared  output  of  the  tide  filter 
was  used.   The  ratios  computed  from  these  integrations  are 
also  given  in  Table  II. 
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From  examination  of  Table  II,  it  was  determined  that 
a  value  of  1/60  for  the  fading  time  constant  gave  the  most 
desirable  results  when  the  integration  was  performed  on  the 
squared  output  of  the  tide  filter.   From  inspection  of  the 
output  statistics  of  the  first  differences,  the  standard 
deviation  of  the  first  differenced  pressure  data  was  ap- 
proximately eight.   Hence,  a  three  sigma  value  of  25  was 
chosen  as  the  threshold  value  of  the  spurious  data  gate 
in  the  first  differencing  subroutine.   Figure  4  shows  the 
effect  of  the  spurious  data  gate  on  the  detector  output 
when  this  threshold  limit  was  set  at  25.   Figure  4  also 
represents  a  fading  time  constant  value  of  1/60  applied 
to  the  integration  of  the  squared  output  of  the  tide  fil- 
ter.  (Compare  Figures  2  and  4.) 

The  critical  threshold  value  of  the  unusual  event  de- 
tector was  determined  to  be  three  times  that  of  the  back- 
ground signal  of  Figure  4(e).   Hence,  the  three  events 
detected  in  Figure  4(e)  will  cause  the  recorder  to  be  ac- 
tivated. 

The  power  spectra  [4]  of  the  tide  filter  and  detector 
outputs  are  shown  on  Figures  5  and  6,  respectively.   In 
Figure  5  the  null  power  readings  coincide  with  even  cph , 
and  the  peaks  coincide  with  odd  cph.    This  result  agrees 
with  the  characteristics  of  the  tide  filter.   The  nodes  of 
the  power  spectral  plot  of  the  detector  output  do  not  ex- 
actly agree  with  those  of  the  tide  filter  because  of  the 
effect  of  the  integration  process. 
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FIGURE  4.  (a)  Pressure,  (b)  First  Differenced  Pressure  (Spurious 
Data  Gate  On),  (c)  Tide  Filtered  Pressure,  (d)  Squared 
Tide  Filtered  Pressure,  (e)  Integrated  Pressure  Versus 
Time    in   Hours . 
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